%% VFA T1 Measurement Simulator with Noise
% Author: Mathieu Boudreau
% Date Created: March 3rd 2012
% Date Modified: April 17th 2012 by Mathieu Boudreau
%

%% Code overview 
% Note: This code in it's current form uniquely implemements the
% linear-least squared procedure of fitting VFA data.
%
% It also requires B1 maps from the AFI method.

%% Clear MatLab 
%

clear all
close all
clc

%%%
tic
%%%

%% Import files
%

% VFA DATA
rawVFAdata(:,:,:,1) = loadminc('REG_SEQ1.mnc');
rawVFAdata(:,:,:,2) = loadminc('REG_SEQ2.mnc');
rawVFAdata(:,:,:,3) = loadminc('REG_SEQ3.mnc');
rawVFAdata(:,:,:,4) = loadminc('REG_SEQ4.mnc');
rawVFAdata(:,:,:,5) = loadminc('REG_SEQ5.mnc');

% AFI DATA
rawAFIdata(:,:,:,1) = loadminc('REG_SEQ6.mnc');
rawAFIdata(:,:,:,2) = loadminc('REG_SEQ7.mnc');


%% Fitting flags
% nlsFlag is set here for future possibility of non-linear least squared
% fitting of the AFI data. Currently (April 2012), nlsFlag = 1 is cannot be
% implemented.

% nlsFlag variable
% 0 = don't do NLS fitting, 1 = do NLS fitting
nlsFlag = 0;

%% VFA Sequence and Image parameters 
% All time sequence parameters are in seconds

seqParam.TR = 0.015;
seqParam.FlipAngles = [3,10,20,30,40];
seqParam.xPixelDim = 2;
seqParam.yPixelDim = 2;
seqParam.zPixelDim = 5;

acqVFAParam.xDim = length(rawVFAdata(:,1,1,1));
acqVFAParam.yDim = length(rawVFAdata(1,:,1,1));
acqVFAParam.zDim = length(rawVFAdata(1,1,:,1));
acqVFAParam.numAngles = length(rawVFAdata(1,1,1,:));

%% Image sizing variables
%

xFOV = [1:acqVFAParam.xDim]*seqParam.xPixelDim;
yFOV = [1:acqVFAParam.yDim]*seqParam.yPixelDim;
zFOV = [1:acqVFAParam.zDim]*seqParam.zPixelDim;

%% Set the AFI input sequence parameters.
% Note that TR1 = N*TR2.

afi.TRratio=5;          % TRratio
afi.NominalAngle=60;    % Nominal flip angle (degrees)

%% Display simulation parameters
%

fprintf('VFA Parameters\n');
fprintf('\n'); %Empty line
fprintf('TR = %1.3f s\n', seqParam.TR);
fprintf('Flip Angle: %1.2f %1.2f %1.2f %1.2f %1.2f\n',seqParam.FlipAngles)
fprintf('\n'); %Empty line

fprintf('AFI Parameters\n');
fprintf('\n'); %Empty line
fprintf('TR ratio = %1.3f\n', afi.TRratio);
fprintf('Nominal Flip Angle = %1.3f\n', afi.NominalAngle);
fprintf('\n'); %Empty line

%% Calculate B1 Map
%

b1Map=AFIB1Map(rawAFIdata(:,:,:,1),rawAFIdata(:,:,:,2),afi.TRratio,afi.NominalAngle);

imagesc(yFOV, xFOV, b1Map(:,:,length(b1Map(1,1,:))/2));
axis image;
title('3.0 T AFI B1 map')
colorbar

%% Calculate Flip Angle Maps
%

flipAngleMaps = zeros(acqVFAParam.xDim, acqVFAParam.yDim, acqVFAParam.zDim , acqVFAParam.numAngles);
for k = 1:acqVFAParam.numAngles
    flipAngleMaps(:,:,:,k) = seqParam.FlipAngles(k);
    flipAngleMaps(:,:,:,k) = flipAngleMaps(:,:,:,k).*b1Map;
end

%% T1 calculation from LLS procedure 
% Note that T1 values smaller than 0 s or larger than 4 s are dumped to 
% T1 = 0s.

[llsT1Map, llsEquiliMagnMap, xVFAFitTerm, yVFAFitTerm] = fitLLSVFAdata(rawVFAdata, flipAngleMaps, acqVFAParam, seqParam.TR);
llsT1Map(llsT1Map>4)=0;
llsT1Map(llsT1Map<0)=0;
llsT1Map=real(llsT1Map);

%% T1 maps
% 

figure()
imagesc(yFOV, xFOV, llsT1Map(:,:,length(llsT1Map(1,1,:))/2));
axis image
title('3T VFA T1 map - zCenter')
colorbar
xlabel('yDirection')
ylabel('xDirection')


figure()
imagesc(xFOV, zFOV, rot90(squeeze(llsT1Map(:,length(llsT1Map(1,:,1))/2,:))));
axis image
title('3T VFA T1 map - yCenter')
colorbar
xlabel('xDirection')
ylabel('zDirection')

figure()
imagesc(yFOV, zFOV, rot90(squeeze(llsT1Map(length(llsT1Map(:,1,1))/2,:,:))));
axis image
title('3T VFA T1 map - xCenter')
colorbar
xlabel('yDirection')
ylabel('zDirection')


%% Plot Whole-volume T1 histogram
%

reshapedllsT1Map = reshape(llsT1Map,1,(acqVFAParam.xDim*acqVFAParam.yDim*acqVFAParam.zDim));
reshapedllsT1Map(reshapedllsT1Map==0)=[];
figure()
[y x]=hist(reshapedllsT1Map,500);
plot(x,y./(length(reshapedllsT1Map)))
title('Whole-volume T1 histogram')
xlabel('T1 (s)')


    
%% End of procedure
%


toc
%%%